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Abstract 

In this article, we design Asymptotic-Preserving Particle-In-Cell methods for the Vlasov-Maxwell 
system in the quasi-neutral limit, this limit being characterized by a Debye length negligible compared 
to the space scale of the problem. These methods are consistent discretizations of the Vlasov-Maxwell 
system which, in the quasi-neutral limit, remain stable and are consistent with a quasi-neutral model 
(in this quasi-neutral model, the electric field is computed by means of a generalized Ohm law). The 
derivation of Asymptotic-Preserving methods is not straightforward since the quasi-neutral model is a 
singular limit of the Vlasov-Maxwell model. The key step is a reformulation of the Vlasov-Maxwell 
system which unifies the two models in a single set of equations with a smooth transition from one to 
another. As demonstrated in various and demanding numerical simulations, the Asymptotic-Preserving 
methods are able to treat efficiently both quasi-neutral plasmas and non-neutral plasmas, making them 
particularly well suited for complex problems involving dense plasmas with localized non-neutral regions. 

Keywords: Plasma, Debye length. Quasi-neutrality, Vlasov-Maxwell, Asymptotic-Preserving scheme. 

1 Introduction 

In a plasma, the Coulomb interaction between charged particles tends to restore the charge neutrality, 
while the thermal motion tends to disturb it. These opposing phenomena introduce a typical length of 
the separation between the electron (ue) and the ion (n^) densities, called the Debye length, and a typical 
oscillation period of the electrons, called the (electron) plasma period. These parameters depend essentially 
on the density and the thermal velocity of the particles. When the scales of interest are large compared to 
the Debye length, the charge separations may be neglected. In other words, the plasma may be assumed 
quasi-neutral. In that case, the Poisson equation is meaningless for the computation of the electric fielcG 

^For the quasi-neutral regime investigated, this also applies to the Maxwell-Ampere equation. 
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E, because this field varies on much larger scales than the charge density. It is thus preferable to assume 
the charge neutrality (this is the so-called plasma approximation) and to compute the electric field by other 
means. Examples of such quasi-neutral descriptions can be found, for instance, in the Langmuir model 
[261138] , the electric field being given by the Boltzmann approximation, or in the Magneto-Hydro-Dynamic 
(MHD) models [7] and the kinetic quasi-neutral models described in [?T1[551 [T|. for which the electric field is 
provided by an Ohm law. As outlined in m chapter 3], “m a plasma it is possible to assume Hi = Ue and 
V • E ^ 0 at the same time” , leading to the following guideline: “Do not use Poisson’s equation to obtain E 
unless it is unavoidable !” which defines the path followed in the present work. 

In many plasma problems, the non-neutral areas are very localized but crucial for the global evolution 
of the plasma. Furthermore, the extent of these regions evolves with time. That is the case, for instance, of 
Plasma Opening Switches (POS) [SHI [SH [JT] where an electromagnetic wave interacts with a dense plasma, 
leading to the formation of non-neutral sheaths. Quasi-neutral models are not appropriate to describe all the 
complex physics that occurs in these non-neutral regions. Non-neutral models, such as the Vlasov-Maxwell 
model or the Euler-Maxwell model, are not appropriate either because the quasi-neutral limit is a singular 
limit (some equations degenerate in the quasi-neutral limit). In particular, explicit discretizations of non¬ 
neutral models are subject to stability conditions on the mesh size and the time step that are all the more 
restrictive than the Debye length and the plasma period are small. 

A first way to cope with these multi-physics problems is to split the domain into non-neutral and quasi¬ 
neutral areas, and to use suitable models in each area (see [15] for example). Another way is to develop 
an Asymptotic-Preserving (AP) method. Introduced by S. Jin |34| for transport in diffusive regimes, the 
Asymptotic-Preserving methods are specifically designed for singular perturbation problems. For a general 
presentation of the Asymptotic-Preserving methodology and examples of application, we refer to |21j . In 
the case of plasma descriptions and their quasi-neutral limit, an Asymptotic-Preserving discretization should 
satisfy the following three properties. 

PI. It is consistent with a non-neutral model. 

P2. It remains stable in the quasi-neutral limit. 

P3. It is consistent with a quasi-neutral model in the quasi-neutral limit. 

These properties make the Asymptotic-Preserving methods able to treat efficiently both quasi-neutral and 
non-neutral areas of a plasma. When the mesh size resolves the Debye length, the Asymptotic-Preserving 
methods provide standard discretizations of a non-neutral model; conversely when the Debye length is not 
resolved, they afford discretizations of a quasi-neutral model. Asymptotic-Preserving discretizations have 
already been developed for various non-neutral models: the Euler-Poisson m, Euler-Maxwell [53], Vlasov- 
Poisson [giia and BGK-Vlasov-Poisson [19] models. The purpose of the present article is to carry on the 
work initiated in these former realizations and to derive Asymptotic-Preserving Particle-In-Cell methods for 
kinetic description of magnetized plasmas, hence the Vlasov-Maxwell model. 

The first step in the design of the Asymptotic-Preserving schemes consists in specifying the quasi-neutral 
model. With this aim, the Vlasov-Maxwell equations are scaled so that they depend on a single dimensionless 
parameter A quantifying how close to quasi-neutrality the plasma is. This parameter A represents different 
diemensionless parameters including the ratio of the Debye length to the space scale of interest and the 
ratio of the plasma period to the time scale of interest. The scaling assumptions are actually similar to the 
assumptions used to derive the MHD models, the kinetic description and the finite inertia of the particles 
being however preserved. In the quasi-neutral model obtained with this scaling, the electric field is computed 
by means of a generalized Ohm law. 

The quasi-neutral model is a singular limit of the Vlasov-Maxwell model. Consequently, a straightforward 
discretization of the Vlasov-Maxwell model would not be consistent with the quasi-neutral model in the 
quasi-neutral limit (condition P3). The key idea is to reformulate the Vlasov-Maxwell system in a set of 
equivalent equations for which the quasi-neutral model is a regular limit. The reformulated system unifies 
the two models, with a smooth transition from one to another according to the value of the asymptotic 
parameter A. By discretizing the reformulated equations, we would obtain schemes that satisfy automatically 
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conditions PI and P3. However, it is preferable to start from the standard Vlasov-Maxell equations, which 
ensures the condition PI, and to use the reformulated equations as a guideline to obtain schemes that 
meet also the condition P3. This approach allows us to use discretizations with good and well-known 
properties. The reformulated equations show that the source term of the Maxwell-Ampere equation must 
be predicted using a discretization of the generalized Ohm law in which the electric field is made implicit. 
We consider two different discretizations of the generalized Ohm law: the first one, namely the AP-Moment 
discretization, is an Eulerian approximation; the other, the AP-Particle discretization, relies on a partial 
Lagrangian approximation (using an advance of the particles). Coupled with a suitable particle pusher, 
these discretizations prove to be stable in the quasi-neutral limit (condition P2). 

Particle-In-Cell (PIC) methods are well documented to lack consistency with the Gauss law, which can be 
the source of non-physical results in the numerical simulations HIMlIl]. There are several ways of enforcing 
the Gauss law in Particle-In-Cell methods (see [3] for a thorough review). In the present work, we focus on 
the elliptic correction of the electric field, which is a simple and robust procedure, and we adapt it to the 
Asymptotic-Preserving framework. 

The Asymptotic-Preserving methods are semi-implicit and share some similarities with the Direct Implicit 
[371 [HI [H [32] and Implicit Moment methods [H [ini [Ml 113 SSI- However, the aim of these methods and 
the way they are derived are quite different. These methods are designed to be free from the usual stability 
constraints on the explicit methods, in order to study large-scale phenomena, but without the aim to ensure 
the consistency with a well-identified asymptotic model. 

The organization of the paper is the following. In Section [31 we scale the Vlasov-Maxwell system, then 
identify the quasi-neutral model and finally reformulate the Vlasov-Maxwell in a set of equations for which the 
quasi-neutral model is a regular limit. These operations are performed first on the standard Vlasov-Maxwell 
system, then on the Vlasov-Maxwell system with elliptic correction. Section [31 is devoted to the derivation 
of the AP-Moment and AP-Particle schemes. Reference explicit schemes are also presented to illustrate the 
defects of the standard discretizations of the Vlasov-Maxwell equations. In Section l4Tl without the purpose 
to be exhaustive, this topic being the subject of an active research for many years, the AP schemes are 
compared with other semi-implicit or implicit PIC methods for the Vlasov-Maxwell equations. In Section 
14.21 the AP schemes are compared with the AP schemes developed for the Vlasov-Poisson equations in [22] . 
Finally, in Section]^ the AP schemes are tested on various and demanding simulations: the classical Landau 
damping; the expansion of a plasma slab into vacuum; a one-dimensional model of POS; the propagation of 
a KMC wave in a two-dimensional model of POS. 


2 The Vlasov-Maxwell system and its quasi-neutral limit 

2.1 The Vlasov-Mcixwell system 

For simplicity, the ions are supposed to form a motionless and uniform background density, denoted by 
rii- The electron evolution is described using a distribution function / depending on the space variable 
X G fix C the microscopic velocity v G ily C and the time t G M"*". The electron density n, 
the electrical charge and current densities, p and J, as well as the stress tensor S are defined from the 
distribution function by 

n= f{x, v,t) dv , p = e{ni — n), J = —e / f{x, v, t)v dv , ^ ^vdv, 

e denoting the elementary charge. The distribution function / satisfies the Vlasov equation 

dtf + v-Vxf--iE + vxB)-V^f = 0, (1) 

m 
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where m is the electron mass, E the electric field, and B the magnetic field. The electric and magnetic fields 
are created by the particles (self-consistent fields) and satisfy the following Maxwell equations: 


^dtE- 

- V X B = —poJ , 

(2) 

dtB + \ 

7 X A = 0, 

(3) 

< 

b] 

II 

P_ 

(4) 


eo 

<1 

ta 

II 

0, 

(5) 


where c is the speed of light, fj,o the vacuum permeability and eg the vacuum permittivity. The above 
equations O-® form the so-called Vlasov-Maxwell system. Of course, this system must be supplemented 
with initial and boundary conditions specific to each problem. Several examples of problems with their initial 
and boundary conditions are presented in Section devoted to numerical simulations. 

In the Vlasov-Maxwell system, the Maxwell-Gauss equation ® (or Gauss law) and the Maxwell-Thomson 
equation ([S]) are actually consequences of the other three equations. The integration of the Vlasov equation 
® over riy gives the continuity equation 

dtp + V-J = 0, (6) 

which translates the conservation of the electron and ion densities. Gombining this continuity equation with 
the divergence of the Maxwell-Ampere equation ([2|), we obtain 

dt(v-E-= 0. (7) 

Taking the divergence of the Maxwell-Faraday equation ®, we find 

dtV ■B = 0. (8) 

The above two equations show that the Gauss law and Maxwell-Thomson equation hold true for t > 0 as 
soon as they hold true at t = 0. 


2.2 Scaling of the Vlasov-Maxwell system 


Gonsider a plasma characterized by the following dimensional parameters: xq for the space scale, to for 
the time scale, no (taken equal to Ui) for the density, vo (taken equal to xo/to) for the velocity, Vthp for 
the electron thermal velocity and {Eo,Bo) for the electromagnetic field. The Debye length Ad (the typical 
length of charge separation) and the (electron) plasma period Tp (the typical period of electron oscillations) 
are defined by: 


Xd = 


e^n-o 


Ty = 


mep 

e^no 


As explained in the introduction, a plasma is considered to be quasi-neutral when the Debye length is very 
small compared to the space scale at which the plasma is observed. Thus, the parameter 



Xo 


quantifies how close to quasi-neutrality the plasma is. 

To investigate the quasi-neutral limit (A —?► 0), we scale the Vlasov-Maxwell system so that it depends 
only on the parameter A. Using the dimensionless variables 


x* = —, t* = -, v* = —, = n* = —, J* 

Xo to Vo no/Vo no 

F* = — B* = — 

Eo' Bo' 


J 

enpvo ’ 
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the Vlasov-Maxwell system can be rewritten as (dropping the stars for the sake of readability): 


dtf + j^v ■ V,/ - vM{E +^vxB)-VJ = 0, 
xB) = -aV, 

P^dtB + V X £; = 0 , 

X^r]M'^y ■£ = 1-n, 

V B = 0, 


(9) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 


where 


M = 


vq 

Vthfi 


7 ] = 


exoEo 


Vo 


= -, P = 


Vo Bo 


En 


The parameter M is the Mach number, 77 the ratio of the electric energy to the drift energy, a the ratio 
of the typical velocity to the speed of light and /3 the induction electric field to the typical electric field. The 
scaling relations defining the quasi-neutral regime are very similar to the most common assumptions of MHD 
models: A <C 1, a 1, M = 1, 77 = 1 , and /? = 1 . A vanishing dimensionless Debye length A —>■ 0 provides 
the quasi-neutrality assumption thanks to the Gauss law. The choice M = 77 = 1 means that the drift 
energy, the thermal energy and the electric energy remain the same order of magnitude. Furthermore, this 
choice implies Tp/to = A, which means that the parameter A controls the smallness of the plasma period with 
respect to the time scale of the problem. Therefore, this quasi-neutral regime is a low frequency asymptotic, 
which explains the stability issues on the time step encountered by standard explicit discretizations. The 
identity ^ = 1 is related to the so-called “frozen-field” assumption, translating the property that, in a dense 
plasma, the magnetic field is convected with the plasma flow. The quasi-neutral regime is also inconsistent 
with the propagation of electromagnetic waves at the speed of light, hence a <C 1. This defines another 
small scale, which is identified to A in the sequel, so that the re-scaled system becomes: 


dtf + vVJ-{E + vxB)-yj = 0, 

(14) 

X^dtE-V X B = -J, 

(15) 

dtB + V xE = 0, 

(16) 

X^V ■E = l-n, 

(17) 

<1 

to 

II 

0 

(18) 


2.3 The quasi-neutral model 

The scaling of the Vlasov-Maxwell system allows us to identify (formally) the quasi-neutral limit of the 
Vlasov-Maxwell system. Setting A = 0 in (fHll - lT^ . we obtain 


dtf + -{E + vxB)-y,f = 0, 

(19) 

V xB = J, 

(20) 

dtB + V xE = 0, 

(21) 

n = 1, 

(22) 

<1 

to 

II 

0 

(23) 


The singular nature of the quasi-neutral limit appears clearly: both the Maxwell-Ampere equation (12011 
and the Gauss law (12^ degenerate. The electric field E can no longer be computed explicitly with the 
degenerate Maxwell-Ampere equation (l20l) . Taking the time derivative of (l20l) together with the curl of 
(EU, we obtain the identity 

VxVxE = -dtJ, (24) 
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which provides a means of computing the solenoidal part of the electric field. But the Gauss law, which is 
degenerate, does not allow us to compute the irrotational part of the electric field. 

While the evolution of the electric field is governed by the displacement current in the Maxwell regime 
(i.e. when the sources do not dominate), it is given by the particles current in the quasi neutral limit. In 
this latter regime, an auxiliary equation is required to derive an Ohm law and explicit the dependence of the 
right hand side of (1241) with respect to the electric field. It is thus necessary to use moments of the Vlasov 
equation (as in [Ml [H [231 Elisa [30]). Multiplying the Vlasov equation by v, then integrating over we 
find the following relation on the current density, which can be viewed as a generalized Ohm law. 


dtJ = ■ S + nE - J X B . (25) 

Using this relation together with (1321) . we obtain an equation that determines entirely E: 

E + yxyxE = JxB-V-S. (26) 

Finally, provided that V x S = J at the initial time, the quasi-neutral model can be (formally) rewritten 
as 

dtf + v-y,f-{E + vxB)-V,f = 0, (27) 

E + VxVxE = JxB-V-S, (28) 

dtB + V X E = 0, (29) 

n = 1, (30) 

V-B = 0. (31) 


Remark 2.1. The rigorous derivation of the quasi-neutral limit for the Vlasov-Maxwell system, i.e. the 
convergence of the solutions of the Vlasov-Maxwell system to a solution of the quasi-neutral Vlasov- 

Maxwell system (HU-dMl) when A —^ 0, is an open problem. Even for the simpler Vlasov-Poisson system, 
few results exist; see for instance the introduction of iSOf for a review. 

2.4 Reformulation of the Vlasov-Maxwell system 

The purpose of this section is to manufacture a set of equations, equivalent to the original Vlasov-Maxwell 
system (I14p - (I18I) . for which the quasi-neutral limit is regular. In other words, the quasi-neutral model (12711 - 
(1311) must be recovered when A is set to zero in this new set of equations. To obtain such equations, we 
reproduce on (fT4ll - (IT^ the operations performed in the previous section on (fT^ - (l23l) . 


The time derivative of (TTKI) together with the curl of (flbll yield 

X^dfE + V X V X U = -dtJ . (32) 

Then, combining the above relation with the generalized Ohm law ()25|) , we find 

X^d^E+ nE + V X (y X E) = J X B-V ■ S . (33) 

Finally, the reformulated Vlasov-Maxwell system is 

dtf + v-V,f-iE + vxB)-Vy = 0 (34) 

X^d^E + nE + V X y X E) = J X B -V ■ S, (35) 

dtB + V X E = 0, (36) 

Xy ■E = l-n, (37) 

V • R = 0. (38) 


This system is equivalent to the Vlasov-Maxwell system (ll4D - (fT51) provided that the Maxwell-Ampere equa¬ 
tion (fTSl) is satisfied at the initial time. 
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2.5 Enforcement of the Ganss law 


In standard Particle-In-Cell methods, the source terms of the Maxwell equations, namely the discrete charge 
and current densities, are computed using a particle-to-grid assignment. The grid quantities do not satisfy 
the discrete equivalent of the continuity equation and thus the Gauss law is not enforced, which can be the 
source of non-physical results in the numerical simulations, as pointed out in [SJi. Different procedures 
are successfully used to correct this deficiency; we refer to [3] for a thorough review. Two main approaches 
can be identified. The first one consists in computing a correction of the electric field, this correction 
being the solution of an elliptic equation (or, in some variants, a parabolic or hyperbolic equation). The 
second approach modifies the particle-to-grid assignement so that the charge and current densities satisfy 
the discrete continuity equation. It is interesting to focus on the first approach, and more specihcally on the 
elliptic correction, because the Gauss law, used to compute the correction, degenerates in the quasi-neutral 
limit. Therefore an adaptation of the elliptic corection needs to be provided. 

A rigorous way to add the elliptic correction is to consider a generalized formulation of the Maxwell 
equations in which the Gauss law is explicitly enforced by means of a correction held (a Lagrange multiplier) 
[1]. The generalized formulation that we use here, also called Boris correction [5], is slightly different, though 
equivalent, to the one proposed in [3]. This formulation involves an electric held E which does not satisfy 
the Gauss law, a corrected electric held E which does satisfy the Gauss law and a correction held p: 


5t/ + u • V,/ - (A + u X H) • V,/ = 0 , 

(39) 

X^dtE-V xB = -J, 

(40) 

dtB + V xE = Q, 

(41) 

X^X7 -E = 1-n, 

(42) 

<1 

to 

II 

o 

(43) 

E = E — Vp. 

(44) 


The boundary conditions prescribed on p are chosen to be compatible with those on E. The generalized 
formulation is well-posed even if the continuity equation is not satished at each time t and is equivalent to 
the standard Maxwell system as soon as the continuity equation is satished at each time t. Gombining (1421) 
and (|44)) , we deduce that p is the solution of the elliptic equation 

\^/\p = X^V ■E-{l-n). (45) 

When the continuity equation holds true at each time t, the generalized equations (I5^ - (H11) are equivalent to 
the standard equations (IT4l) - (fT8l) . Indeed, in this case, the electric held E satishes the Gauss law A^V • E = 
1 — n. Therefore, p vanishes and E = E. 

At the time-discretized level, the generalized formulation provides us with an obvious way to compute an 
electric held satisfying the Gauss law. First, an electric held is computed using (14(1)1 - (|4T|) : then, this electric 
held is corrected using and (H51) . 

If we add the elliptic correction to the quasi-neutral model, we obtain 

E + V X {V X E) = J X B - V ■ S , 
dtB + V X E = Q, 
n = 1 , 

V • B = 0, 

E = E — Vp. 

The degeneracy of the Gauss law prevents from computing p directly. To obtain an equation on p, we need 
again to use the moments of the Vlasov equation. However, the continuity equation is assumed to be not 
exactly satished by the moments of the distribution function. Following the spirit of the Boris correction. 
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this unconsistency is corrected in the equation by means of an electrostatic deviation of the electric field, 
giving rise to the following modified continuity equation 


dtJ = V ■ S+ n{E-Vp) - J X B , (46) 

dtn = y-J, (47) 

Introducing these definitions into the time derivative of the Gauss law yields to 

-X^dfAp - V • {nVp) = -E-W^ :S-W- {nE) + V • (J x S), 

which is finally equivalent, assuming that the Gauss law is satisfied at initial times, to 

- X^d'^Ap - V • (nVp) = : S -V ■ (nE) + V • (J x S). (48) 

Remark 2.2. The source term of this equation is the difference beetween the time evolution of the density 
produced by the Vlasov equation and that of the electric field divergence predicted thanks to the moments. 
This point will be clarified with the time semi-discretization introduced in the sequel (see section \3.3.1\) . The 
equation (14811 provides a means of preserving the consistency with the Gauss law whatever the values of X 


and is therefore AP in the quasi-neutral limit. 

The quasi-neutral model with correction can thus be stated as follows 

dtf + v-V^f-{E + vxB)-V,f = 0, (49) 

E + VxVxE = JxB-V-S, (50) 

dtB + V xE = Q, (51) 

Ap = -d(n + -.S + V -E-V -{J xB), (52) 

V-R = 0, (53) 

E = E-Vp. (54) 

Similarly, the reformulated Vlasov-Maxwell system with correction can be rewritten as 

dtf Ev-V^f -{E + vxB)-V,f = Q, (55) 

X^d(E + nE + VxVxE = JxB-V-S, (56) 

dtB + V X E = Q, (57) 

X^d(Ap + V • (nVp) = -d(n -b : 5 -f V • [nE) - V • (J x R), (58) 

V-R = 0, (59) 

E = E-Vp. (60) 


3 Asymptotic-Preserving schemes 

Now that the Vlasov-Maxwell system has been scaled and the quasi-neutral model has been identified, we 
can state rigorously the properties that an Asymptotic-Preserving discretization must satisfy. 

PI. For A > 0, the discretization is consistent with the Vlasov-Maxwell system (l39l) - ((4^ (or (|55]) - (l60l) l. 

P2. The stability conditions on the time step and the mesh size do not depend on A. 

P3. For A = 0, the discretization is consistent with the quasi-neutral model (H^ - (I5^ . 




In this section, we derive two Asymptotic-Preserving schemes, called AP-Moment scheme and AP-Particle 
scheme. They use a Yee finite-difference approximation for the fields |50j (with a regular rectilinear grid) 
and standard assignment-interpolation procedures for the particles, such as the nearest grid point or cloud- 
in-cell procedures 0133]. Other choices could have been made: finite volumes or finite elements instead 
of finite differences for instance. The schemes are presented in a three-dimensional spatial setting and, for 
simplicity, the domain is assumed to be a rectangular parallelepiped with periodic boundary conditions. It 
is straightforward to derive one-dimensional and two-dimensional versions of the schemes. 

We also present two reference explicit schemes. These schemes illustrate the defects of the standard 
discretizations of the Vlasov-Maxwell equations and will be compared to the Asymptotic-Preserving schemes 
in the numerical tests. Finally, we derive 

3.1 Definitions and notation 

3.1.1 Discrete fields and discrete vector calculus operators 

We consider different kinds of discrete fields on the grid (which is rectilinear and regular): primal and dual 
scalar fields, edge scalar field, primal and dual vector fields, primal symmetric second-order tensor field. 

• The values of a primal scalar field are located at the vertices of the cells, while the values of a dual 
field are located at the centers of the cells. The values of an edge scalar field are located at the center 
of the edges. 

• The components of a primal vector field are located at the center of the edges: the x-, y-, and z- 
components are located at the edges oriented in the a;-, j/-, and z-direction, respectively. The compo¬ 
nents of a dual vector field are located at the center of the faces: the x-, y-, and z-components are 
located at the faces normal to the x-, y-, and z-direction, respectively. 

• The diagonal components of a primal symmetric second-order tensor field are located at the vertices 
of the grid. The xy-, xz- and yz-components are located at the center of the faces normal to the 
x-direction, y-direction, and z-direction, respectively. 

Discrete differential operators can be defined on the discrete fields defined above by using central finite 
differences (and assuming periodic boundary conditions). 

• A discrete curl operator V;iX is defined for the primal and dual vector fields. When applied to a primal 
vector field (resp. dual vector field), the discrete curl operator yields a dual vector field (resp. primal 
vector field). Furthermore, if Fh is a primal vector field and Gh a dual vector field, then 

VhXFh-Gh-VhXGh-Fh=0. (61) 

• A discrete divergence operator V/i- is defined for the primal and dual vector fields. When applied to a 
primal vector field (resp. a dual vector field), the discrete divergence operator yields a primal scalar 
(resp. a dual scalar field). If Fh is a primal or a dual field, then 

V/, • (Vh X F,,) = 0 . (62) 

• A discrete gradient operator V/j is defined for the primal and dual scalar fields. When applied to a 
primal scalar field (resp. a dual scalar field), the discrete gradient operator yields a primal vector (resp. 
a dual vector field). 

• A discrete divergence operator ^/h- is defined for the primal tensor field. It yields a primal vector field. 

• A discrete cross product x h between a primal vector field and a dual vector field is defined. It yields 
a primal vector field. Such a discrete operator is built using local averages. 
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3.1.2 Notation 


• The grid spacings in the a;-, y-, and z-direction are denoted by Ax, Ay, and Az, respectively. Let 

h = ^/ 

• The time interval is discretized with a uniform time step At. Let V = 7 At, for any 7 £ 1R+. 

• The discrete electric and magnetic fields at time t'’' are denoted by E'^ and The discrete electric 
field is a primal vector field while the discrete magnetic field is dual vector field. The discrete correction 
field, denoted by is a primal scalar field. 

• Let N be the number of particles. The vectors containing the position and the velocity of the particles 
at time E are denoted by and V^, respectively. The position and velocity vectors of the jth 
particle at time E are denoted by ^ and j, respectively. 

• The value of a field Fh interpolated at the position Xnj is denoted by Fh{XNj). 

• The discrete electron density accumulated from the particles at position Aat as a primal scalar field 
(resp. an edge scalar field) is denoted by nhiX^) (resp. hh{XN)). The discrete current accumulated 
from the particles of position Xn and velocity Vn as a primal vector field is denoted by Jh{XM, Vn)- 
The discrete second-order moment accumulated from the particles of position Xjq and velocity Vn as 
a primal tensor is denoted by Sh{XN, Vn)- 


3.2 Reference explicit schemes 

The first reference scheme combines a Boris scheme for advancing the particles and a leap-frog discretization 
for the Maxwell equations (SJ Chapter 15]: 


where 


At " ’ 

Vje{i,.. 

-.N}, 

(63) 

Y+ W” 

^ N ,o ^ ^ N ,j jjm f vm 

2 X Bf^ (Aa,j), 

Vje{i,.. 

■,N}, 

(64) 

^m+1 _ T^m 1 1 

\2^h V7 sz 7 ( vm-\-l 

^ At ,Vj^ ), 



(65) 

^ X Tir = 0 , 



( 66 ) 

^ ^m+1 _ ^ 



(67) 

X^AhPh = - (1 - nhix^+^)) , 



( 68 ) 


Vje {i,...,A}, 


(69) 

yN,j = VnP + \^tE^ , 

Vje {i,...,A}, 


(70) 




(71) 


The above equations are solved in the order (IMll . (l64l) . (l63l) . (1651) . (l68ll . (l67ll . so that the scheme is fully 
explicit. This scheme is subject to a number of stability conditions that are given below (between brackets, 
the conditions are expressed with dimensional parameters and without the scaling assumptions). The time 
step must resolve the plasma period: 


At < 2A [At < 2Tp] . 


(72) 
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The grid spacing must resolve the Debye length: 


h<(X [h< CXd] , 


(73) 


where C is a parameter depending on the assignment-interpolation procedure. Otherwise, aliasing will heat 
up the plasma. Furthermore, the time step and the grid spacing must satisfy a Courant condition involving 
the velocity of the electrons. 


T<^ 


At ■ 

vthfi-r < 1 
h 


(74) 


and and another one involving the speed of light. 




(75) 


The first three constraints are due to the explicit discretization of the particle motion, while the last constraint 
is due to the explicit discretization of the Maxwell equations. Note that the constant in the right-hand side of 
the Courant condition ((75ll is specific to the Yee finite-difference discretization [50] . Because of the stability 
conditions ([72]) and (ESI), the above scheme is unstable in the quasi-neutral limit. Therefore, it does not 
satisfy the condition P2, and hence is not Asymptotic-Preserving. 

In the second reference scheme, the leap-frog discretization of the Maxwell equations (l65|) - (l66l) is replaced 
by an implicit 0-scheme (with ^ < 0 < 1) but the sources remain explicit: 


,E\ 


'm +1 


-EV 


At 

r^m+l T^rn 

At 


- V. X B 


m-\-6 _ T / \rm-\-l 


V. X = 0 


(76) 

(77) 


with 

= eE^+^ -h (1 - e)E]^ , = 0B™+1 + {1- 0)B^ . 

The 0-scheme is unconditionally stable for ^ < 0 < 1. Therefore, this second scheme is not subject to 
the stability condition H75D . However, it is still subject to the conditions (17^ and (1751) . and hence is not 
Asymptotic-Preserving. The computation of and requires the solution of a linear system, which 

makes this scheme slightly more costly than the first one. Other implicit schemes than the 0-scheme could 
have been used; see for instance mum- The properties of the 0-scheme (stability, energy conservation, 
dispersion) are recalled inlAl 


3.3 Asymptotic-Preserving schemes 

3.3.1 General framework 

To derive AP schemes, it is preferable to start from the standard Vlasov-Maxell equations (l3^ - (l44l) than 
from the reformulated equations (1551) - (l60l) . It allows us to use discretizations with good and well-known 
properties. The reformulation operations described in Sections 12.41 and 12.51 are however used as a guideline 
to obtain schemes consistent with the quasi-neutral model. In this section, we derive the common general 
structure of two AP schemes, called AP-Moment and AP-Particle schemes. Their specificities are addressed 
in the next two sections. 

1. First, the Maxwell equations are discretized using a 0-scheme (with j < 0 < 1): 

pm+1 _ rpm 

\2^h V7 s/ f>En-\-6 _ jm+1 

^ - - ~'Jh 

pTTi+l _ TDm 

I V7 N/ Tj'm-\-6 _ n 

-—- +VhXEj^ - 0 , 


(78) 

(79) 
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where 


E]^+<^ = + (1 - e)E]^ , B]^+^ = 9B'^+^ + (1 - 9)B^ . 

An implicit discretization of the Maxwell equations is needed to avoid a Courant condition similar to GSI. 
The current used as source term in the discrete Maxwell-Ampere equation (175)) is defined by an approx¬ 
imation of the generalized Ohm law (1251) . It is crucial to make the electric field implicit in this approximation 
to ensure the consistency with the quasi-neutral model. We consider two types of approximation. The first 
one, called AP-Moment, is based on an Eulerian integration. The second one, called AP-Particle, relies on a 
partial Lagrangian approximation (using an advance of the particles). In both approximations, the current 
can be written in the form 


rm+l _ 7-m+l,* 

Ju — Ou 


■AtnhiX^)E', 


m\ pm+1 


h 


(80) 


where is a quantity depending on the scheme. Note that and are primal vector fields. 

The equations ((78l) - (l80l) are equivalent to the linear system 




Urv 




where A™ is a linear operator defined by 


■^h 


Eh 

Bh 


VAt2 


nh{X^))Eh-^yhXBh 


At 2 


Bh 


At 


XhX Eh 


(81) 


(82) 


for any fields Eh and Bh- In practice, the electric field E^~^^ and the magnetic field B^'^^ are computed 
by solving this linear system. It is invertible and well-conditionned even when A <C At and A <C h. Indeed, 
using the identity (1611) . we obtain the inequality 



> 


i— 

I At 2 


min {nh{X^)) ] \\Eh\\ 


At 2 


WBhf 


When A <C At and A <C h, the density nh{X'^) is expected to be close to 1 (the charge separation is negligible 
at scales far larger than A). Therefore, the quantity [}?/At^) -|- min (nh{X^)) is close to 1. 

2. In a second step, the electric field is corrected. Equations (IT5]) and (Hill are discretized as follows: 

X^Xh ■ = 1 - , (83) 

E^+l ^ ^m+l _ 

The charge density used as source term in (1831) is built using an approximation of the moment equations 

(H51) and (1371): 


= nh{X^) + AtXh ■ Jh^^, (85) 

J-+1 = J-+i’* + Ath^(A]0)£;;(*+i. (86) 


In the above equations, is a primal scalar field and is a primal vector field. It is essential to make 

the electric field implicit in this approximation to guarantee the consistency with the quasi-neutral model. 
From ((551) - (1551) . we deduce that the correction ph is solution of the equation 


-V 


h ’ 


At 2 


+ nh{XN) ^hPh = 


^ l-nhiX^) _^^JfX 


At 2 


V At2 


■MX. 


N. 


E 


m+1 


-AtXh-Jh^^'*. (87) 


12 







Just as the linear system (ISTl) . the above linear system (157)) is invertible and well-conditionned even when 
A <C At and X h. Using the divergence of (l78l) . it can be rewritten in a simpler form; 


- V 


h ' 


A^2 + '‘V Njj^hPhj ^^2 ^ 


( 88 ) 


Note that the right-hand side of (l88)) evaluates the inconsistency of the Gauss law at time t™ which amounts 
to the difference beetween the density accumulated from the particles and the electric field implicitly predicted 
at the previous time step. In practice, the correction field ph is computed with (1881) . then the electric field 
is corrected with (l84l) . 

3. Finally, the particles are advanced with a Boris-like scheme: 


Vje A}, (89) 

-U-+1 ^ X Vj e {1,..., A}, (90) 

where 

= yN,V - , vj e {1, ..., A}, (91) 

yN,j = yN,j + l^tEi:^^H^N,j), vje{i,...,A}. (92) 

With this particle pusher, the scheme is subject to the Courant condition (1741) . However it is not subject 
to conditions dUl) and dzi and thus satisfies Property P2. Moreover, it presents favorable conservation 
properties. The choice of the particle pusher is discussed in more detail iniBl 


^N,j 

At 

T/m+l T^m 

At 


3.3.2 AP-Moment scheme 


The AP-Moment Scheme relies on the following approximation of the generalized Ohm law (l25l) : 


im+l 

Ot 


At 


-Xh-Sh{X^,V^) 


h^(A; 0 )A™+i + U^) = 0 . 


(93) 


This amounts to set 

= j^(a;;),ix) + At(v^ • Sh{x's, u^) - mx's, v^) x,, . (94) 

Proposition 3.1. The AP-Moment Scheme, defined by (1751) - (1501) . (15^ . (1551) . (1501) - (1051) and (10^ . is con¬ 
sistent with the Vlasov-Maxwell system when A > 0. It is consistent with the quasi-neutral model (I49l) - (l54p 
when A = 0. 


Proof. For simplicity, we fix 0 = 1. First, the system (15T1) - (1551) yields 


At2 V 


£,^+1 _Elf] = 




-X/hXX/hX - hh{X^)Ejf+^ - Vffi) X Bff 


The Ampere law is assumed to be satisfied at the precedent time step, with 

A 2 


V, X sr - vffi) ^ ^ - K~") , 


so that the following approximation holds 


— f 

At2 V 


^m +1 _ 


Vft, X Vh X E' 


m +1 

h 




'm+1 I V 7 cm jm ( vm -\rm\ .. tditi 
h +yh- S^ - Jf, (Ajv, ) X 


( 95 ) 
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defining a consistent discretization of the reformulated Ampere law (1341) . 

Next, the equations (l83]) - (l86l) and (1^ provide 

■ ^ (l - V, ■ Ar) 

+ • J^X^, v;^) + V^ 5 "* - X B^) - Vh ■ {nh{X'S)K^^) ■ 

Assuming that the Gauss law as well as the continuity equations are satisfied at the previous time step so 
that the following approximations hold 

X^Vh-E^^l-nh{X^). 

AtVn ■ JJTiX^. V^) ^ , 


the following equation can be stated 

- ^^ ■ ((^ + MX^))VhPH) ^ ^ + 2 n,(A-) - nu{X^-^)) 

- Vl : 5"* - V, • (jr X B^) - V, • , 

which defines a time discretization of the reformulated Gauss law (l58)) provided that the correction at time 
level m and m — 1 vanish. 

□ 


3.3.3 AP-Particle scheme 

In the AP-Particle scheme, the contributions oi J x B and V • S' in the generalized Ohm law (1^51) are 
approximated with a particle advance. Observing that 


V • S- J X S 



{v ■ Xxf - {vx B) ■ V„/) V dv , 


(96) 


we set 


tTTI-I- 1, X / 1 T x777.-l-1 


(97) 


with 


^N,j ^N,j 

At 

T/m 

''n,j 

At 


■\r7Tl -\-1 




Vj e {1,. 

..,iV}, 

(98) 

Vj G {1,. 

..,iV}. 

(99) 


Compared to the AP-Moment scheme, the AP-Particle scheme requires an additional push of the particles, 
from the state {X^,V^) to the state {X^^^'* However, it avoids the computation of the tensor 


Sh. 


Proposition 3.2. The AP-Particle Scheme, defined by (1751) - (15(11) . (15^ . ((551) . (1551) - (1M1) and (1571) - (1551) . is 
consistent with the Vlasov-Maxwell system when A > 0. It is consistent with the quasi-neutral model (HSD-dMl) 
when A = 0. 


Proof. When the number of particles is sufficiently large, the following appoximation is valid: 
J^(A;;;+l’^ V^+i’*) A(A™, V^) -AtVn- Sh{X^,Vffi) + AtJh{X^,Vffi) X 


( 100 ) 


This proves that the AP-Particle scheme shares the same consistency properties than the AP-Moment scheme. 

□ 
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4 Comparison with other implicit PIC methods 

4.1 Fully implicit, Direct Implicit and Implicit Moment methods 

Some level of implicitness is required in the discretization of the Vlasov-Maxwell equations to obtain AP 
schemes. It is needed to ensure the consistency with the quasi-neutral model and the stability in the 
quasi-neutral limit. Other implicit kinetic methods have been developed since the 1980s. Designed for the 
simulation of large-scale phenomena, they are derived to relax the main stability conditions that explicit 
methods must satisfy. We refer to [33] for a recent review. In this section, we discuss the AP character of 
the three main classes of implicit methods (the fully implicit methods, the Implicit Moment methods and 
the Direct Implicit methods) and compare them with the AP-Moment and AP-Particle schemes. 

Theoretically, a fully implicit discretization of the Vlasov-Maxwell system is stable for any discretization 
parameters. The resulting problem is a huge system of coupled nonlinear equations (the particle equations 
and the field equations). Impressive realizations have been achieved in the past few years with the use of 
Jacobian-Free-Newton-Krylov solvers, preconditioning, optimized and often massively parallel implementa¬ 
tions Nevertheless, fully implicit discretizations are still too costly for multi-dimensional 

simulations. In fully implicit discretizations, the electric field at the advanced time level occurs in the 
discretization of the terms dE jdt and V x A in the Maxwell equations and is also involved in the definition 
of the source terms, via the particle equations. Therefore the consistency with the quasi-neutral model is 
recovered when the coupled system is solved. 

The Implicit Moment methods |4l[Tni|53lll5l|5lll6l|40lll3] are semi-implicit methods. They decouple 
the held equations from the particle equations by using macroscopic evolution equations on p and J to 
predict the sources of the held equations at the advanced time level and This approach reduces 

dramatically the computational cost compared to fully implicit methods, while retaining favorable stability 
properties. The use of moment equations to predict the sources at the advanced time level is a common point 
with the AP-Moment and AP-Particle schemes. In particular, the discretization of these moment equations 
is very similar in the Implicit Moment methods and the AP-Moment scheme. 

The Direct Implicit methods [n ETjlHlIsg are also semi-implicit. They enjoy the same accuracy and 
stability properties as the Implicit Moment methods. They rely on a linearization of the fully implicit 
discretization, which decouple the field equations from the particle equations. The sources of the field 
equations are linearized around explicitly extrapolated positions of the particles. The source prediction 
in the Direct Implicit methods can actually be interpreted as an approximation of the moment equations, 
comparable to the one used in the AP-Particle scheme. The extrapolation step plays essentially the same role 
as the first particle advance in the AP-Particle scheme (though being more intricate in the Direct Implicit 
approach). 

From this brief review, we can infer that the fully implicit. Implicit Moment and Direct Implicit methods 
are generally AP in the quasi-neutral limit. The fully implicit methods are far more costly than the AP- 
Moment and AP-Particle schemes. The Implicit Moment and Direct Implicit methods share some similarities 
with the AP-Moment and AP-Particle schemes in their formulation. However, their motivation and the 
methodology used to derive them differ significantly. The aim of the AP-Moment and AP-Particle schemes is 
to be consistent with a clearly defined quasi-neutral model, not to relax stability conditions. The elimination 
of the stability conditions related to the plasma period and the speed of light is a consequence of scaling 
assumptions made in Section 12.21 for the definition of the quasi-neutral limit (the scaled plasma period and 
the ratio of the typical velocity to the speed of light vanish in the quasi-neutral limit). The derivation of the 
AP-Moment and AP-Particle schemes relies on a reformulation of the Vlasov-Maxwell system which unifies 
the Vlasov-Maxwell model and the quasi-neutral model in a single set of equations. This reformulation 
highlights the terms that need to be built or implicited in order to ensure consistency with the quasi-neutral 
model. This methodology allows us to limit the computational cost of the schemes to what is necessary in 
view of the AP property. Furthermore, this methodology could be applied to quasi-neutral models with a 
more reduced complexity, leading to more efficient numerical methods for some kinds of problems. 
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4.2 AP-Moment and AP-Particle schemes in the electrostatic regime 

The electrostatic regime is characterized by a vanishing magnetic field. In the dimensionless system (ED-dlSl), 
this amounts to the asymptotic /3 —?► 0. In this regime, the electric field is irrotational, since the Maxwell- 
Faraday equation m simplifies into V x FI = 0, and it is assumed to derive from a scalar potential (there 
exists a scalar field (p such that E = —Vp). As a consequence, the Gauss law is sufficient to determine 
completely the electric field and the Vlasov-Maxwell system reduces to the Vlasov-Poisson system: 


5t/ + u-V,/ + V<^-V„/ = 0, (101) 

-A2A(/)=l-n. (102) 


Let us examine the AP-Moment and AP-Particle schemes in the electrostatic regime. Setting the magnetic 
field to zero, they read 


+ hh{X'S)]E, 


7^m+l 


-T 

- X 


1 

^ rpm _ rm + l,* 

^ At ^ 


+ hh{XN))XhPh) = 


1 - n^(A-) A- 


At 2 


-A^v.-Ar, 


"N,j 


At 


N,j _ Trm+1 
~ ^N,j : 


VjG {1,...,A}, 


•trm+l J/m 

At 


= -E 


r+'(x]o,^), vjg{i,...,a}, 


with, for the AP-Moment scheme. 


= J^(X™, VPP) + AtVh- Sh{X^, VJP). 


and, for the AP-Particle scheme, 


j7Tl-\- \ .'k 


Mxp^ + Atvpp,v;p). 


(103) 

(104) 

(105) 

(106) 

(107) 

(108) 

(109) 


In a one-dimensional spatial setting, the AP-Moment and AP-Particle schemes can be simplified further. 
Indeed, for any primal vector field Eh, there is a primal scalar field ph such that Eh = —Vph- Let and 

be the potentials associated with Ejp~^^ and respectively. Equation (I105|) implies the equality 

Xhp'h^^ = Xh +Ph^- Then, adding (11031) and (I104p . we deduce 


This equation determines entirely the potential and thus the electric field Ejp~^^. Therefore, the 

Maxwell-Ampere equation can be disregarded and the schemes reduce to the equations (11061) . (11071) and 
(IllOp . We remark that the AP-Moment scheme is, in this case, equivalent to the PIC-AP2 scheme introduced 
in [22]. 


5 Numerical simulations 

In order to assess their efficiency and investigate their properties, the AP schemes are tested on various 
problems and compared with reference explicit schemes (those described in Section 13.2|) . The first two 
problems, namely the classical Landau damping and the expansion of a plasma slab into vacuum, are one¬ 
dimensional in space and velocity and purely electrostatic. The two other problems, which describe a POS 
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(Plasma Opening Switch), involve magnetized plasmas. The first one is a very simplified model of POS (one¬ 
dimensional in space and two-dimensional in velocity); the second one, more realistic, is two-dimensional in 
space and velocity and shows the propagation of a KMC wave. 

The framework of the above problems is not exactly the general framework described in Sections [5] and 
[3] (the problems are one-dimensional or two-dimensional, the magnetic field is sometimes disregarded, the 
ions are not always motionless). The numerical schemes are implemented accordingly. These developments 
being straightforward, we do not detail them. However, for sake of clarity, we state the version of the 
Vlasov-Maxwell system nsed for each problem. All the simulations are performed with the cloud-in-cell 
assignement-interpolation method and the parameter 9 is always taken equal to 1 in the 0-schemes. 

5.1 Landau damping 

When a plasma in a spatially homogeneous equilibrium state is slightly perturbed, it returns exponentially 
fast, with oscillations, toward its initial equilibrium state. This is the so-called Landau damping. This 
problem allows us to test the ability of the AP schemes to accurately reproduce phenomena occurring at the 
Debye length and plasma period scales. 

The non-neutral model used for this problem is a one-species, one-dimensional, electrostatic model (the 
magnetic field is disregarded, the ions are motionless and form a uniform background). Scaled as in Section 
the equations are 


dtf + vdxf - Edvf = 0, (111) 

X^dtE = -J, X'^dxE = l-n. ( 112 ) 

The space domain is [0,47r] and the scaled Debye length A is taken equal to 1. The initial electron density 
follows a Maxwellian distribution with a small spatial perturbation: 

/o(a;,u) = (l-tacos (D)(113) 

where a = 5 • 10“^. Periodic and homogeneous Dirichlet boundary conditions are prescribed for the particles 
and the electric field, respectively. An analytical approximation of the electric field can be computed by 
applying a Laplace-Fourier transform to the linearized system; see [20] for the detailed calculation. Keeping 
only the dominating mode, the others being quickly damped, the following approximation is obtained: 

Eix,t) « -1.4708Q;e-°-i®33*cos(1.4156t- 0.536245) sin (a;/2). (114) 

Numerical simulations are performed with discretization parameters smaller than the Debye length and 
the plasma period. The numerical results, together with the analytical approximations given by the formula 
(I114|l . are represented in Figured) We observe that the AP schemes reproduce quite accurately the oscil¬ 
lation period of the Landau damping but they are more dissipative than the reference explicit scheme (the 
first one described in Section Note that the discrepancy between the numerical results and the ana¬ 
lytical approximation at the first oscillation is due to the neglected Laplace-Fourier modes in the analytical 
approximation. The conservation properties of the AP schemes for the plasma oscillations could probably be 
improved with a more careful design of the time discretization, with centered or high-order approximations 
instead of backward approximations. 

5.2 Plasma expansion 

The second problem is the expansion of a plasma slab IMlIllllI]- The plasma is initially confined in a 
small area at the center of the domain and surrounded by vacuum. A non-neutral sheath forms at the 
plasma-vacuum transition and the large electric field created in this sheath accelerates the ions, leading to 
the plasma expansion. While the quasi-neutral model is able to account for the ion motion and the plasma 
expansion in the plasma bulk, it is inaccurate in the sheath. Therefore, this problem allows us to verify. 
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Figure 1: Landau damping. Analytical approximations (Ana.) and results computed using the reference 
scheme (Ref.), the AP-Particle scheme (AP-P), the AP-Moment scheme (AP-M). The simulations are per¬ 
formed with Ax ~ A/20, At = A/10 and 10® particles. 


on the one hand, the consistency of the AP schemes with the quasi-neutral model and their stability in the 
quasi-neutral regime and, the other hand, the consistency with the non-neutral plasma description, with the 
transition from one regime to the other. It is also an excellent test for the energy conservation properties 
of the schemes since it relies on a kinetic energy transfer from the electrons to the ions via the electric field 
created in the sheath. 

The non-neutral model used for this problem is a two-species, one-dimensional, electrostatic model (the 
magnetic field is disregarded, the ions are not motionless). This test case is stated and implemented in 
dimensional variables, the equations being: 

dtfe + V d^fe - —E dyfe = 0, 

nie 

dtfi + V dxfi + —E dvfi = 0, 
rrii 

dtE = --, dxE=^, 
eo eo 

where the indices i and e denote the quantities related to the ions and the electrons, respectively. 

The space domain is [—L,L] with L = 1 m. The problem being symmetric, the simulations are actually 
performed only on the half domain [0,L]. The ion mass is such that rrii = 1836 me. The initial ion density 
is equal to riio for x € [0,0.02] and is zero in the rest of the domain (different values of riio are used in the 
simulations). The initial electron density Ueo satisfies the Maxwell-Boltzmann relation Ueo = n-io exp((/)o), 
where t/o is the electrostatic potential, solution of — A^o = e(nio — neo)/eo (this nonlinear problem is 
solved numerically). The initial electron and ion velocities follow Maxwellian distributions with zero mean 
velocities and respective temperatures Tg and T^. The electron temperature Tg is chosen so that the initial 
electron thermal velocity is equal to Vth,e = 1 m-s“^ and the ion temperature is such that Ti = 10“®Te. To 
represent the symmetry at the left end of the half domain, a specular reflection condition is prescribed for 
the distribution functions (i.e. exiting particles are reinjected with reversed velocities) and an homogeneous 
Dirichlet condition is enforced on the electric field. At the right end, an absorbing condition is enforced on 
the distribution functions (i.e. exiting particles are not reinjected) and an homegeneous Neumann condition 
is enforced on the electric field. In addition to the time and space scales related to the electron motion, this 
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problem involves time scales related to the ion motion, namely the ion plasma period Tpi = 

and the ion accoustic wave speed Cs = The speed of the plasma expansion is theoretically the 

same order as [55]. 




(a) Normalized ion density ni/nio at time t = SOr^i, as a (b) Normalized electric field Ef y/{nokBTe,o)/eo at time 

function of the space variable x/\]j. t = 30rpj, as a function of the space variable xIXd. 



Ref. Wt 
Ref. We 
••■^■Ref. Wi 
-^Ref. We 
AP-P Wt 
-G AP-P We 
-e-AP-P W, 
-e-AP-P We 
AP-M Wt 
AP-M VPe 
-V-AP-M Wi 
-<^AP-M We 


(c) Total (Wt)-, electron kinetic (We)-, ion kinetic (Wi) and electric (We) energies (normalized by the initial total energy) as 
functions of the time ijr-pi. 


Figure 2: Low-density plasma expansion. Results computed using the reference scheme (Ref.), the AP- 
Particle scheme (AP-P) and the AP-Moment scheme (AP-M). The simulations are performed with Aa; = 
0.4 Ad, At = O.bTp and 10® particles for each species. 


A first series of simulations is performed with a low-density plasma (the initial density riio is adjusted to 
obtain Ad = 10 “® m and Tp = 10 “® s). The mesh size and the time step are smaller than the Debye length and 
the plasma period. The numerical results are represented in Figure|2]and, to facilitate the comparisons, they 
are expressed with the same normalized quantities as in [28l [22]. The AP schemes reproduce correctly the 
physics of the plasma expansion: the plasma reaches the position xt ~ 136Ad at the end of the simulation, 
there is a large electric held at the plasma edge, the (thermal) electron kinetic energy is converted into (drift) 
ion kinetic during the simulation. The total energy of the system is conserved, which demonstrates the good 
conservation properties of the schemes. Furthermore, the results of the AP schemes are almost identical to 
the results of the reference explicit scheme (the hrst one described in Section [3T2ll and are in agreement with 
the results of [551155] . 

A second series of simulations is performed with a high-density plasma (the initial density riio is adjusted 
to get Ad = 10“^ m and Tp = 10“^ s). The reference explicit scheme is still used with discretization 
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parameters resolving the Debye lenth and the plasma period, whereas the AP schemes are now used with 
discretization parameters significantly larger than the Debye length and the plasma period (the number of 
particles is also reduced). The Maxwell-Boltzmann relation used in the previous simulations to compute 
the initial electron density turns out to be inaccurate with a coarse mesh. It is thus better to make the 
initial electron density equal to the ion density. The numerical results are represented in Figure El As 
in the previous simulations, the AP schemes reproduce correctly the physics of the plasma expansion and 
the energy conservation is satisfactory. The results provided by the AP schemes are even quite close to 
the results obtained with the reference explicit scheme, despite the important difference of computational 
cost. Note that, in Figure |3(a)[ the oscillations in the reference scheme results are not physical and are 
due to the insufficient number of particles. The outputs of the AP-Moment scheme and the AP-Particle 
schemes are similar but the discretization used for the AP-Moment scheme is finer. This shows that the 
AP-Particle scheme is less diffusive than the AP-Moment scheme and provides a better tracking of the 
plasma-vacuum interface. Finally, these simulations demonstrate the ability of the AP schemes to simulate 
quasi-neutral problems with moderate computational costs. However, it is not possible to use too coarse 
discretizations, owing to the deterioration of the energy conservation. This is a well-known drawback of the 
semi-implicit methods |25) . due to fast particles crossing more than one cell in a time step. This point should 
be investigated in subsequent realizations. 

5.3 A one-dimensional model of POS 

A Plasma Opening Switch (POS) is a device used to deliver a large current with a rapid increase of its 
impedance |24) . It consists of a coaxial cylindrical transmission line filled whith a high-density plasma and 
connected to an input power generator [55]. In a first phase, the conduction phase, the plasma short-circuits 
the two electrodes of the transmission line and prevents the power to be delivered to the load. Then, the 
interaction of the electromagnetic wave with the plasma leads to the formation of a vacuum gap and finally 
to the total opening of the plasma, making possible the transmission of the current to the load. During 
the opening of the plasma, charge separation phenomena are essential and non-neutral sheath appear at the 
plasma edge 1221147]. 

A very simplified model of POS, one-dimensional in space and two-dimensional in velocity, is considered 
in this section. This model has been introduced in [23]. Let us denote by {x,y,z) the three-dimensional 
Cartesian coordinates. The particles move only along the x-direction but have a velocity both in the x- 
direction and the j/-direction so that the electromagnetic field generated by the particles has only components 
Ex, Ey and B^- Therefore, the two-species Vlasov-Maxwell model reduces to the following equations: 

dtfe + Vx ■ dxfe - —{Ex + VyBz)dv^fe - —{Ey- VxBz)dv fe = 0 , 

rrie rrie 

dtfi + Vx ■ dxfi + —{Ex + VyB^)dy^fi + —{Ey - VxBx)dy fi = 0, 

rrii mi 

-^dtEx = -yoJx , -^9tEy + dxB,, = -yoJy , dtB,, + dxEy = 0, 

supplemented with the Maxwell-Gauss and Maxwell-Thomson equations. In the above equations, the indices 
i and e denote the quantities related to the ions and the electrons, respectively. 

The space domain is [0, L] with L = 0.2 m. At the initial time, the plasma fills the region between 
X = 0.05 and x = 0.15 with the density no (different values of no are used in the simulations). The initial 
ion and electron velocities follow Maxwellian distributions with the same temperature Ti = Te = 10^ eV. 
The ion mass is such that = 2 • 10^ mg. A transverse electromagnetic wave is sent from the left end of 
the domain. The value of its electric component at x = 0 and time t is 

Ey,inc{t) = Ainc {t/tinct^^ (s + {t/Uncf ((V^inc)® “ 

where Ainc = 1-8 ■ 10® V-m“^ and fine = 10“® s. Transparent boundary conditions are prescribed at each 
end of the domain to avoid wave reflections. 
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at time t = 233rpi, as a function of the space variable 
x/Xd- 



(b) Ion mean velocity (normalized by Cs) at time t = 
233rpi, as a function of the space variable xjXj:). 



(c) Total (Wt), electron kinetic (We), ion kinetic (Wi) and electric (We) energies (normalized by the initial total energy) as 
functions of the time 


Figure 3: High-density plasma expansion. Results computed using the reference scheme (Ref.), the AP- 
Particle scheme (AP-P) and the AP-Moment scheme (AP-M). The discretization parameters are : Ax = 
0.5 A_d, At = 0.2Tp and lO”^ particles (for each species) for the reference scheme ; Ax = 25A£), At = 
and 10® particles for the AP-Particle scheme ; Ax = IOAd, At = 2rp and 10® particles for the AP-Moment 
scheme. 


Four series of simulations are performed (their main characteristics are collected in Table [Ij. In the Low-a 
simulations, the AP schemes and the reference explicit scheme (the second one described in Section [SA]) are 
used with discretization parameters that resolve the Debye length and the plasma period. The numerical 
results, represented in Figure |4l are comparable for the three schemes. In particular, the level of numerical 
noise is equivalent. We observe that the electrons and the ions are accelerated by the incident wave (see 
Figures |4l(c) and|4l(d)). The electrons being less massive than the ions, they are accelerated more strongly 
and are expelled from the plasma edge, which breaks the quasi-neutrality in this zone (see Figures |4l)a) and 
SKb)). This charge separation creates a large electric field (see FigureHKe)). As for the magnetic field, it 
is significantly transmitted through the plasma (see FigureUKf)). 

The Low-b simulations deal with a higher initial plasma density. The reference explicit scheme still uses 
discretization parameters that resolve the Debye length and the plasma period, while the AP schemes use 
a mesh size 20 times larger than the Debye length and a reduced number of particles. Despite the huge 
difference of computational cost between the simulations, the numerical results are indistinguishable (see 
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(c) Electron current density in the x-direction and y- (d) Ion current density in the x-direction and in¬ 
direction (in A-m“^). direction (in A-m”^). 




Figure 4: Low-density POS (Low-a). Results computed using the reference scheme (Ref.), the AP-Moment 
scheme (AP-M) and the AP-Particle scheme (AP-P) at time t = 3 ■ 10“® s. 
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Table 1: One-dimensional model of POS. Plasma parameters and discretization parameters used in the 
different simulations. The initial density no is in m“^, the initial Debye length A_d in m and the plasma 
period Tp in s. The value Np is the total number of particles. 


Config. 

no 

A_d 


Reference scheme 

AP schemes 

Tp 

Ax/Xd 

<1 

Np 

Ax/Xd 

At/xp 

Np 

Low-a 

IQi® 

10-4 

o 

7 

o 

1—1 

1 

0.1 

10® 

1 

0.1 

10® 

Low-b 

1017 

5 10-® 

5 10-11 

0.2 

0.04 

107 

20 

0.04 

104 

High-a 

102O 

10“® 

10-11 




10^ 

1 

10® 

High-b 

1022 

10-7 

10-1® 




104 

10^ 

10® 


Figure [S]). In these simulations, the current created by the electron motion at the plasma edge is strong 
enough to stop the penetration of the magnetic field into the plasma bulk (see Figure [SKd)). 

For the high-density cases (High-a and High-b), the computational cost of the reference explicit scheme 
is prohibitive, so that no simulation is carried out with this scheme. The AP schemes are used with very 
coarse discretizations. For the most demanding case, the mesh size is 10^ larger than the Debye length, 
the time step is 100 times larger than the plasma period and the number of numerical particles is only 
10®. The numerical results are shown in Figure [6l In contrast to the low-density simulations, the incident 
magnetic field is almost entirely reflected at the plasma edge (due to the very large electron current). The 
level of numerical noise in the outputs remains moderate, even in the High-b simulation. These simulations 
demonstrate the ability of the AP schemes to handle high-density plasmas and vacuum-plasma transitions 
with very coarse discretizations. 


5.4 Propagation of a KMC wave in a POS 

We now consider a two-dimensional model of POS where the plasma density is not uniform in the transverse 
direction to the electromagnetic wave propagation. In this configuration a magnetic shock wave, the so- 
called KMC wave, propagates into the plasma [49]. The existence of KMC waves can be derived from the 
quasi-neutral equations (see El). Consequently, their simulation offers a means to verify the consistency with 
the quasi-neutral limit in an electromagnetic context. 

The model is a one-species, two-dimensional model. The ions are motionless but their density is not 
uniform. Let us denote by (x, y, z) the three-dimensional Cartesian coordinates. The particle motion is re¬ 
stricted to the (x, 2 /)-plane and thus the electromagnetic field generated by the particles has only components 
Ex, Ey, and B^- Finally, the Vlasov-Maxwell system simplifies into 


dtf + Vx ■ dxf + Vy ■ dyf - —{Ex+ VyB^)dyJ - —{Ey - VxB^)dy f = 0 , 

m m 

^dtEx - dyB^ = -fioJx , ^9tEy + dxB^ = -fioJy , -H dxEy - dyEx = 0 , 


(116) 

(117) 


supplemented with the Maxwell-Gauss and Maxwell-Thomson equations. 

The space domain is the rectangle [0,Lx] x [0, Ly] with Lx = 0.2 m and Ly = 0.03 m. The lower and 
upper sides of the domain represent the cathode and the anode of the POS, respectively. At the initial time, 
the plasma fills the area between x = 0.05 and x = 0.15 with the density 


no{x,y) = < 


T^min'^max / [( ^min ^max) ( 2 //O.OI 1) -)-Un 


^ ^min 


if 0 < y < 0.01 
if 0.01 <y < 0.02 
if 0.02 <y< 0.03 
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(b) Electron current density in the a:-direction and y- 
direction (in 




Figure 5: Low-density POS (Low-b). Results computed using the reference scheme (Ref.), the AP-Moment 
scheme (AP-M) and the AP-Particle scheme (AP-P) at time t = 4 • 10“® s. 
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(a) High-a. Magnetic field Bz (in Teslas). 
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(d) High-b. Magnetic field Bz (in Teslas). 


Figure 6: High-density POS (High-a and High-b). Results computed using the AP-Moment scheme (AP-M) 
and the AP-Particle scheme (AP-P) at time t = 4 • 10“® s. 
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Table 2: Propagation of a KMC wave in a POS. Plasma parameters and discretization parameters used in 
the different simulations. The inital densities nmin and rimax are in m“^, the electron temperature Te in eV, 
the initial Debye length Ad in m, the theoretical KMC wave speed in m-s“^. The value Np is the number of 
particles. 


Config. 

^min 

^max 

Te 

Ad 

%KMC 

Grid 

Ax 

Ay 

Np 

(a) 

1019 

102O 

6-104 

1.69- 10"® 

25 - 10® 

100 X 100 

2-10-3 

3-10-4 

10® 

(b) 

1019 

1021 

6-102 

5.35- 10-® 

25 - 10® 

400 X 100 

5 - 10-4 

3-10-4 

4-10® 


A transverse electromagnetic wave is sent from the left end of the domain. Transparent boundary conditions 
are prescribed at each end of the domain to avoid wave reflections. 

The initial plasma density is such that dy {x, y) = 10^(ninax —?T-min)/(?T-max'ramin) for 0.01 < y < 0.02. 
Therefore, according to the theory presented in O the incident wave should propagate into the plasma as a 
KMC wave (between the lines y = 0.01 and y = 0.02). Moreover, the expected speed of this KMC wave is 


Vkmc = 10 ' 


> ^max ^mii 
€I-Iq ^max^min 


(118) 


Two simulations are carried out on this problem, both with the AP-Moment scheme (note that the Gauss 
elliptic correction is not implemented in these simulations). The plasma and discretization parameters are 
specified in Table [2] The evolution of the magnetic field for the simulation (a), reported in Figure [7l shows a 
rapid magnetization of the plasma in the region where the density gradient is located. The electrons emitted 
at the cathode produce a current which prevents a uniform penetration of the magnetic field in the plasma. 


This current is gradually deflected with the propagation of the magnetic field, as depicted in Figures 7(c) 
and |7(d)] We evaluate the speed of the KMC waves in the simulations by marking the position of a certain 
magnetic level set at different times (see Figure [8]). The observed values are in agreement with the theoretical 
values: in the simulation (a), it is 76% of the theoretical value, while in the simulation (b), it reaches 82% 
of the theoretical value. These good numerical results are obtained with discretizations that do not resolve 
the Debye length. In particular, for the discretization (a), the mesh size is 10^ times larger than the Debye 
length. 


6 Conclusion 

We have derived two Asymptotic-Preserving Particle-In-Cell methods for the Vlasov-Maxwell system in 
the quasi-neutral limit. The scaling assumptions made for the definition of the quasi-neutral limit, similar 
to those used to derive the MHD models, yield a kinetic quasi-neutral model where the electric field is 
computed by means of a generalized Ohm law. The Asymptotic-Preserving methods are consistent with 
either the quasi-neutral model or the Vlasov-Maxwell model according to how the discretization parameters 
resolve the plasma parameters, which make them able to simulate complex plasma problems with a reasonable 
computational cost. No rigorous numerical analysis is provided in this article and this should be the subject 
of future work. However, the numerous numerical investigations demonstrate conclusively the efficiency 
of the methods to account for phenomena evolving at the plasma period and Debye length scales, as well 
as quasi-neutral phenomena usually well described by the MHD theory. In particular, they are able to 
cope with vacuum-dense plasma interfaces and the formation of non-neutral sheaths. Although other semi- 
implicit Particle-In-Cell methods (the Direct Implicit and Implicit Moment methods) present the same kind 
of Asymptotic-Preserving properties, the methodology developed in this article is new and provides a rigorous 
framework for the quasi-neutral limit problem. Furthermore, it opens the way for addressing more singular 
asymptotics and deriving more efficient numerical methods for some kinds of problems. 
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(b) Magnetic field Bz (in Teslas) at time te = 3.67 ns. 



(c) Current density Je (in A- m iii decimal log- (d) Current density Je (in A- m in decimal log- 

scale) at time = 2.58 ns. scale) at time te = 3.67 ns. 


Figure 7: Propagation of a KMC wave in a PCS. Numerical results computed using the AP-Moment scheme, 
with discretization (b). 
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(a) Configuration (a). Level sets of the magnetic field at times = 1.75, = 2.16, = 2.56, 

4“' = 3.03, 4“^ = 3.56 and 4“^ = 4.09 ns. 



(b) Configuration (b). Level sets of the magnetic field at times = 1.24, = 1.62, = 2.06, 

= 2.58, = 3.11 and 4*’' = 3.67 ns. 
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Figure 8: Propagation of a KMC wave in a POS. Estimation of the KMC wave speed in the simulations (a) 
and (b). The velocities Vi-j are the estimated velocities on the time interval They are computed 

using the position of the level set = —0.8 Tesla. 
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A Yee finite differences and ^-scheme for the Maxwell equations 


We consider the discretization of the homogeneous d-dimensional Maxwell equations (d =1, 2 or 3) with Yee 
finite differences for the space approximation and a 0-scheme for the time integration. Using the discrete 
operators introduced in Section ITTI the discrete equations are 


with 


1 _ T^rn 

V7 N/ + ^ A 

? -- 

pm+1 _ TDm 

^ ^ + V, X = 0 > 


+ (1 - 


-b (1 - e)B^ 


(119) 

( 120 ) 
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The above scheme is second-order accurate in space. It is first-order accurate in time for 9 ^ \ and second- 
order accurate for 0 = The following energy balance holds true: 

, ( 121 ) 

where £J^ = -I- Therefore, the d-scheme is unconditionally stable for 9 G [ 5 , 1 ]- It is 

dissipative for 9 G] 5 ,1] (all the more dissipative than 9 is large) and energy-conserving for 9 = ^. Unlike the 
leap-frog scheme, the 0-scheme is dispersive even in a one-dimensional setting and with a Courant number 
equal to I. 


B Comparison of some particle pushers for the AP schemes 


We examine the properties of some particle pushers within the AP schemes. The particle pushers we consider 
are variants of the Boris scheme: 


T^m+1 _ ym 

^N,j ^N,j T^m+6 

At ’ 

ym+1 _ 

- iXN,j ) 




W, 


N,j 


X Bl 


I ( \rm \ 


( 122 ) 

(123) 

(124) 

(125) 


with (a, b, c) € { 0 , 1 }^. 

Numerical simulations show that the electric field and the velocity must be made implicit (5 = I, c = 1) 
to overcome the stability condition (|7^ . They also show that, whatever the choice for (a,5, c), the scheme 
remains subject to the Courant condition (1741) . An explicit discretization of the position (a = 0) is preferable 
to an implicit discretization (a = I), since it yields a scheme easier to solve and more accurate. In particular, 
in the case of a constant electric field and a zero magnetic field, the choice a = 0 yields a symplectic scheme, 
unlike the choice a = I. Symplectic time-integration schemes ensure excellent conservation properties and 
an accurate behavior in long-time simulations [29) . 


C KMC waves 


The existence of KMC waves can be derived from the quasi-neutral equations (with motionless ions) presented 
in Section I 2 T 3 I Neglecting the inertia term dtJ and the pressure term V • S' in the generalized Ohm law (l25ll 
and rewriting it with dimensional variables, we obtain the relation 

nE - J X B = 0 

e 

Then, using the Maxwell-Ampere equation and Maxwell-Faraday equations (still with dimensional variables), 
we deduce that the magnetic field evolution is governed only by the Hall term: 


dtB = —V X ((V X B X B)/{jjLoen)) . 

In the two-dimensional setting ()116|) - (I117|1 . if the density does not vary along the a;-axis and does vary along 
the y-axis, the above equation simplifies into a Burgers-like nonlinear hyperbolic equation: 


dtB, = -- 


2e/ro 


do^Bi 
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This equation admits shock wave solutions, with a speed 



They are the so-called KMC waves, named after Kingsep, Mokhov and Chukbar. 


(126) 
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